clear; clc; close all;

% 模型参数
a = 0.1;
b = 0.1;
k = 1.85;
func = @mptm;

% 扫描参数 c 的范围
c_min = -0.2;
c_max = 1.05;
c_step = 0.001;
c_range = c_min : c_step : c_max;
num_c = length(c_range);

% 迭代设置
transient = 100;   % 舍弃暂态迭代次数
N = 5000;           % 用于 LE 计算及采样的有效迭代次数
total_iter = transient + N;

% 初始状态（这里固定为 (0.1, 0.1)）
x0 = 0.1;
q0 = 0.1;

% 用于存储每个 c 下的 Lyapunov 指数
LE1_vals = zeros(num_c,1);
LE2_vals = zeros(num_c,1);

% 用于分岔图：存储每个 c 下采样得到的 x 值
bifurcation_data = [];  % 每行记录 [c, x]

%% 参数扫描
for i = 1:num_c
    c_val = c_range(i);

    % 初始化状态为固定初值
    x = x0;
    q = q0;

    % 先迭代 transient 次，去除暂态
    for t = 1:transient
        [x, q] = func(x, q, a, b, c_val, k);
    end

    sample_points = 100;
    sample_interval = floor(N / sample_points);

    temp_bif = zeros(sample_points, 2); % [c, x]
    for t = 1:N
        [x, q] = func(x, q, a, b, c_val, k);
        % 如果达到采样间隔，则记录当前 x 值
        if mod(t, sample_interval) == 0
            idx_sample = t / sample_interval;
            temp_bif(idx_sample, :) = [c_val, x];
        end
    end
    % 累计保存分岔数据
    bifurcation_data = [bifurcation_data; temp_bif];

    [LE1, LE2] = LEs(func, x0, q0, a, b, c_val, k, N, transient);
    LE1_vals(i) = LE1;
    LE2_vals(i) = LE2;

end

% 子图 1：Lyapunov 指数
subplot(2,1,1);
hold on; box on; grid on;
plot(c_range, LE1_vals, 'r', 'LineWidth',1.2);
plot(c_range, LE2_vals, 'b', 'LineWidth',1.2);
xlabel('c'); ylabel('Lyapunov Exponent');
title('Lyapunov Exponents vs. Parameter c');
legend('LE1','LE2','Location','Best');
yline(0, 'k--', 'LineWidth',1);

% 子图 2：分岔图
subplot(2,1,2);
hold on; box on; grid on;
scatter(bifurcation_data(:,1), bifurcation_data(:,2), 2, 'k', 'filled');
xlabel('c'); ylabel('x'); 
title('Bifurcation Diagram');

